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FOREWARD 


Numerical  models  for  calculating  ocean  sound 
transmission  loss  continue  to  be  of  high  interest,  es¬ 
pecially  those  which  are  useful  for  calculations  in  an 
environment  that  changes  with  range.  The  model  described 
in  this  report  is  a  powerful  tool  for  application  in  such 
environments . 

When  this  report  was  prepared,  the  author 
worked  for  the  Acoustic  Environmental  Support  Detachment 
(AESD)  of  the  Office  of  Naval  Research.  Prior  to  publi¬ 
cation,  the  function  of  AESD  was  incorporated  into  the 
Numerical  Modeling  Division  (Code  320)  of  the  Naval 
Ocean  Research  and  Development  Activity  (NORDA) .  The 
report's  contents  are  of  sufficient  interest  and  general 
applicability  to  warrant  making  them  available  by  pub¬ 
lication  at  this  time  as  a  NORDA  technical  note.  The 
original  preparation  of  the  material  was  sponsored  by 
the  Long  Range  Acoustic  Propagation  Project  of  ONR  (now 
NORDA  Code  600). 


DARRELL,  CAPT. 
ding  Officer,  NORDA 
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ABSTRACT 


The  Parabolic  Equation  Model  is  a  wave- 
acoustics  model  designed  for  the  computation  of  acoustic 
transmission  loss  as  a  function  of  range  and  depth  in 
range  dependent  ocean  environments.  The  elliptic 
reduced  wave  equation  is  approximated  by  a  parabolic 
partial  differential  equation  that  can  be  numerically 
integrated  by  marching  the  solution  forward  in  range. 

The  model  is  primarily  useful  for  predicting  low  fre¬ 
quency  ( <  200Hz)  acoustic  propagation  of  energy  along 
waterborne  paths.  This  report  briefly  describes  the 
physics  and  mathematics  of  the  model  and  documents  a 
computer  program  developed  at  AESD.  Individual  routines 
are  documented  in  an  appendix.  Environmental  input 
routines  must  be  supplied  by  the  user  and  are  not  des¬ 
cribed  in  this  report. 
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THE  AESD  PARABOLIC  EQUATION  MODEL 


I  INTRODUCTION 

The  Parabolic  Equation  (PE)  Model  replaces  the  reduced 
elliptic  wave  equation  with  a  parabolic  partial  differential 
equation  that  can  be  integrated  numerically  using  the 
Tappert-Hardin  split-step  Fourier  algorithm.  This 
parabolic  approximation  is  useful  when  propagation  is 
predominantly  radial  and  backscatter  can  be  neglected. 

The  parabolic  wave  equation  includes  diffraction  and  all 
other  full-wave  effects  as  well  as  range  dependent 
environments.  The  numerical  algorithm  has  exponential 
accuracy  in  depth,  second  order  accuracy  in  range  and 
is  unconditionally  stable.  The  entire  range  and  depth 
dependent  acoustic  field  is  computed  as  the  solution  is 
marched  forward  in  range.  The  model  assumes  a  flat 
pressure  release  ocean  surface  and  a  vanishing  field  at 
the  maximum  depth  of  the  finite  Fourier  transform.  A 
pseudo  radiation  condition  is  introduced  at  the  water- 
bottom  interface  by  smoothly  attenuating  the  field. 

Since  the  error  in  the  parabolic  approximation  increases 
as  angles  increase  from  the  horizontal,  steep  slopes  can 
cause  inaccuracies.  Also,  computing  time  increases  as 
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frequency  increases.  Thus,  the  model  is  primarily  useful 


I 
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for  predicting  low-frequency  "  <  200Hz"  acoustic  propa¬ 

gation  of  energy  along  waterborne  or  shallow  angle  bottom- 
bounce  paths.  In  addition  the  numerical  solution  is  limited 
by  practical  Fourier  transform  sizes  and  computing  time. 
Therefore  the  frequencies  that  can  be  handled  by  the  model 
depend  on  the  input  environment. 


II  THE  PARABOLIC  APPROXIMATION 

In  a  medium  of  constant  density,  cylindircally 
symmetric  about  an  axis  containing  a  point  time-harmonic 
source,  the  acoustic  pressure  P  satisfies  the  reduced 
elliptic  wave  equation 

[v2  +  K02  (r,z)J  P  (r ,  z)  =  0  (1) 

where  K0  is  a  reference  wave  number 

K0  =  (2) 

Co 

(j  is  the  angular  frequency  of  the  source,  CG  is  a  reference 
sound  speed,  and  n  is  the  index  of  refraction 

co 

n  (r ,  z)  =  -  (3) 

c (r,z) 

Writing  the  acoustic  pressure  in  the  form 

P(r,z)  =  Hr,z)  H0(l)(K0r)  (4) 

where  the  Hankel  function  H0 ^ 1 ^  represents  the  primary  radial 
dependence  of  the  field  in  terms  of  an  outward  propagating 
cylindrical  wave,  and  under  the  assumption  that  the  observation 
range,  r,  is  many  wavelengths  from  the  source  (i.e.  K0r>>l) , 
we  may  make  use  of  the  asymptotic  form 


H0(1>  (K0r) 


i(K0r-  n/k) 


(5) 


to  obtain 


tprr  +  2iK0^r  +  i|/2z  +  K02  [n  (r,z)-l]ip  =  0  (6) 

Up  to  this  point  we  have  made  cylindrical  symmetry  and 
far  field  approximations.  We  now  consider  propagation  that 
is  predominantly  radial  and  make  the  additional  "parabolic1' 
approx imat ion 

<prr  <<  2iK0  \pr  (7) 

which  is  equvalent  to  neglecting  the  reflected  field  and  are 
left  with  the  Leontovich-Fock^  parabolic  equation  for  the 
transmitted  field 


=*  i  (A+B)  tjj 

where  A  and  B  are  the  real  operators 

a  =  _L  !L 
2Ko  3z2 

B  -  ^  (n2-l) 


(8) 


(9) 


(10) 
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and 


(There  is  a  constant  of  proportionality  to  be  determined 
for  ip  wnich  is  related  to  the  source  strength  normaliza¬ 
tion.  See  Section  3.2.)  This  approximation  was  first 

used  in  underwater  acoustic  applications  by  Tappert  and 
2 

Hardin. 

McDaniel"*  examined  the  effect  of  this  parabolic  approx¬ 
imation  on  the  propagation  of  normal  modes  in  a  laminar  ocean 
and  has  shown  that  discrete  modes  are  propagated  with  the 
correct  amplitudes  and  mode  shapes,  but  with  errors  in  phase 
and  group  velocities.  This  error  can  cause  substantial 
shifts  in  the  modal  interference  pattern.  In  general  the 
solution  of  the  parabolic  equation  well  approximates  the 
solution  of  the  elliptic  wave  equation  (equation  (1))  over 
a  narrow  region  of  the  spectrum  of  eigenvalues.  The 
eigenmode  with  range  phase  velocity  corresponding  to  the 
reference  sound  speed  is  propagated  correctly.  Hence 
the  parabolic  approximation  can  be  described  as  a  "narrow- 

4 

band"  approximation. 

In  many  applications  where  a  single  group  of  modes  or 
in  the  ray  analogy,  a  single  ray  family,  dominates  the 
propagation,  the  parabolic  approximation  is  entirely 

4 


5 


adequate.  One  simply  selects  the  proper  turning  point 
velocity  as  the  reference  sound  speed.  However,  when  many 
modes  are  propagating  such  a  choice  is  not  possible  and 
the  resultant  errors  may  be  unacceptable. 

As  a  means  of  improving  the  narrowband  parabolic  approx¬ 
imation,  Buchal5  has  suggested  a  mapping  of  the  index  of 
refraction  that  approximately  preserves  the  phase  velocities 
and  turning  points  of  the  normal  modes  with  turning  points 
in  the  water  column.  The  sound  speed  profile  is  perturbed  to 
produce  a  modified  index  of  refraction  such  that  the  phase 
velocities  of  the  parabolic  solution  modes  equal  the  phase 
velocities  of  the  elliptic  solution  modes.  The  WKB  approx¬ 
imation  is  then  used  to  map  the  mode  turning  points.  The 
resultant  mapping 

(z,n)  -*■  ( z , h )  =  |zn1//z,  (zn-1)  (12) 

significantly  improves  the  agreement  between  the  elliptic 
and  parabolic  solutions . 


6 


2 


Ill 


THE  NUMERICAL  ALGORITHM 


A .  Description 

For  a  useful  long-range  low-frequency  acoustic  propa¬ 
gation  model  we  require  the  numerical  solution  of  equation 
8  for  realistic  range  dependent  ocean  environments.  Several 
numerical  methods  are  applicable  to  the  solution  cf  equations 
of  this  form.  Me  have  chosen  to  implement  the  split-step 
Fourier  method  of  Tappert  and  Hardin. ^  This  algorithm  is 
applicable  to  a  particular  class  of  equations  and  has 
several  advantages.  In  addition  to  exponential  accuracy 
in  Z  and  second  order  accuracy  in  r,  it  is  exactly  energy 
conserving,  unconditionally  stable,  and  computationally 
efficient.  The  algorithm  does,  however,  require,  a  uniform 
mesh  and  periodic  boundary  conditions  in  Z  because  of  its 
use  of  the  fast  Fourier  transform.  The  sound  speed  profiles 
must  also  be  filtered  to  avoid  introducing  high  angular 
frequency  components  into  the  acoustic  field.  This  means 
that  discontinuities  must  be  smoothed. 

7 

Following  the  analysis  of  Buchal  and  Tappert,  a 

solution  of  equation  8  can  be  formally  represented  as 

/r+Ar 

U (?)  dr 

r  1 


+  E 


iK  r ,  z ) 


iArU 

e 


+  i  U 


rr 


(r ,  z) 


(13) 
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where  U  -  A+B 


A  and  B  are  defined  by  equations  9  and  10 
and  U,  Ur,  Urr  are  evaluated  at  r  +  Ar . 

Neglecting  terms  of  third  and  higher  orders,  the  numerical 
algorithm  becomes 

iAr (A+B) 

4>(r+Ar ,7.)  -  e  i|)(r,z)  (14) 


The  implementation  of  the  split-step  algorithm  requires 
the  assumption  that  the  operators  A  and  B  commute,  which 
they  do  not.  If  the  operator  is  written  as  follows 


iAr (A+B)  iArA  iArB  iArA 
2  2 
=  e  e  e 


1  3  J  I[a(AB-BA)  -  (AB-BA)B] 

+  i  [B  (AB-BA)  -  (AB-BA)B] 


(15) 


+  0 


f(Ar)4  ) 


then,  again  to  third  order  in  Ar,  the  numerical  algorithm 
takes  the  form 


iArA  iArB  iArA 

T  2 

^(r+Ar,z)  =  e  e  e  ip  (r,z) 


(16) 
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iArB 


The  evaluation  of  e  is  a  straight-forward  multipli- 

i  ArA 

cation.  Since  A  is  a  differential  operator  e  is 

evaluated  by  means  of  a  fast  Fourier  transform  operation. 
Thus  the  form  of  the  split-step  algorithm  implemented  here 
becomes 


ip  (r+Ar,z) 


iArKn (n2-l) 


FFT 


-1 


-iArK2 

2K0  , 

e  FFT  l\p  (r , 


(17) 


where  FFT  is  the  spatial  Fourier  transform  from  z  to  K  and 
FFT--*-  is  the  inverse  transform.  This  algorithm  advances 
the  field  by  alternating  between  two  stages.  The  first 
stage  advances  the  field  as  if  the  propagation  were  in 
a  homogeneous  medium  thus  accounting  for  the  effects  of 
diffraction.  The  second  stage  then  accounts  for  the 
effects  of  the  environment  on  the  propagation. 

B .  Implementation 

The  essence  of  the  parabolic  equation  technique  is  the 
numerical  solution  of  the  Leontovich-Fock  parabolic  wave 

equation  by  the  Tappert-Hardin  split  step  Fourier  algorithm, 
subject  to  appropriate  initial  and  boundary  conditions.  To 
satisfy  the  pressure  release  boundary  condition  at  the  ocean 
surface 


ij>(r,o)  =  0 


(18) 
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\p(r,z )  is  made  anti-symmetric  in  z.  This  is  accomplished 
by  replacing  the  complex  Fourier  transform  operations  in 
equation  17  by  real  sine  transform  operations  and  transforming 
the  real  and  imaginary  parts  of  the  field  separately.  A 
pressure  release  surface  at  the  maximum  depth  sample 
considered  in  the  transform  is,  however,  implicit  in  this 
procedure.  This  non-physical  boundary  condition  is  handled 
by  the  following  artifice:  We  extend  the  depth  sampling  in 
the  transform  1/4  the  maximum  physical  depth  relevant  to 
the  problem  and  assume  an  index  of  refraction  in  this  region 
of  the  form 

”  /2~zmax\ 

n2  =  n2  +  iae  \  6  )  (19) 

where  nB  is  the  index  of  refraction  at  the  water-bottom 
interface  and  zmax  is  the  maximum  depth  sample  in  the 
transform.  This  form  attenuates  the  additional  discrete 
modes  introduced  by  truncating  the  transform  and,  with 
proper  choice  of  a  and  0,  models  a  pseudo  radiation  condition 
at  the  water  bottom  interface.  The  constants  a  and  6  were 
chosen  empirically  to  give  agreement  with  a  normal  mode  pro¬ 
gram  for  the  ocean  index  of  refraction  over  a  homogeneous  half¬ 
space.  The  values  found  by  this  process  are 
a  =  .01 

0  =  1/3  the  artificial  bottom  depth. 
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Before  the  split-step  algorithm  can  be  used  to  march 


the  field  forward,  the  field  must  be  defined  over  the 
entire  depth  mesh  as  an  initial  condition  at  some  range. 

In  addition,  the  numerical  method  requires  that  this  field 
be  band  limited.  In  the  implementation  considered  here, 
the  initial  field  is  assumed  to  be  a  Gaussian  beam 


-<z-Z0)2 
\l>  (0,  z)  =  S  e 


(20) 


where  S  is  an  effective  source  level,  W  is  a  beamwidth,  and 
Z0  is  the  depth  of  the  point  source  beinq  modeled.  To 
evaluate  S  and  W  we  solve  the  parabolic  wave  equation  in 
a  homogeneous  medium 


2iK0^r  +  1^22  —  0 


(21) 


with  equation  20  as  an  initial  condition  obtaining 


~1/2-  (z-Zo) 


\p  (r , z)  =  S 


l+“£  2 

K0WJ 


WZ/l+i2r  \ 

V  w2) 


(22) 


From  equation  11 


11 


while  from  the  solution  of  the  elliptic  wave  equation  for 
a  point  source  in  a  homogenous  medium 


2  _ 


1+  (Z-ZQ) 


(24) 


Matching  the  parabolic  and  elliptic  wave  equation  results 
for  large  r  yields 


for  the  effective  source  level  of  the  beam  and 


(25) 


(23) 
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2 


(26) 
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for  the  width. 

The  parabolic  phase  velocity  correction  (equation  12) 
is  implemented  as  follows : 

1.  A  uniform  mesh  is  established  in  z  the  transformed 
depth  variable . 

2.  The  sound  speed  is  defined  to  be  piecewise  linear 

in  z 

C(z)  =  Cj+  g  (z-zj) 

3.  The  mapping 


/  2 


implies 


g^j2+  Zj  V~( gZj) 


2  +  4C0  (Cx-gZi) 
2C0 


and 


C 


o 


C,  +  g(Zj-Zi) 


where  CQ  is  the  reference  sound  speed 


4.  The  transformed  index  of  refraction  on  the  uniform 


mesh  is  given  by 

rij  =  (2n j-1)  >/* 

A  1-2-1  filter  is  then  applied  to  the  uniform  mesh  (z,  n) 

(or  (2  ,  n)  if  no  parabolic  phase  velocity  correction  is  to 
be  applied)  to  smooth  the  gradient  discontinuities  associated 
with  the  boundaries  in  the  piecewise  linear  sound  speed 
profile . 

Two  forms  of  error  control  are  employed  during  the  march. 
The  range  step  is  constrained  by  requiring  the  ratios  of  the 
error  terms  (third  order)  in  equation  13  and  equation  15  to 
the  first  order  terms  be  small.  To  evaluate  these  terms  we 
approximate  the  reduced  pressure  locally  as  a  plane  wave 

iKQ  [r  (cos0-l)  +  zSinOl 

<J>(r ,  z)  ~e  (27) 


Equation  14  then  yields  for  the  range  integration  error  term 


Ar<<24 


(n2-l-Sir20)2  +  +  ii%ll£ 


K, 


K0(n2)rr  +  (n2)zzr  +  i  Kq  (n2)rz  Sin9 


(28) 


For  n2  independent  of  r  this  constraint  does  not  limit  Ar 
and  has  not  been  implemented  in  the  current  version  of  the 
computer  code. 


Similarly  from  equation  15  we  have  for  the  first  order 
term 


iAr  (A+E )  =  i^£Ar  (r2-l-Sin2®)ii/ 

2 


(29) 


and  for  the  third  order  term 


i  (Ar) 


V4  [a(AB-BA)- (AR-BA)^f  1  /2  j^B  (AB-BA)  -  (AB-BA)P.J 


(30) 


1K0 

“48 


(Ar) 


((n2)z)2+V- 


(n) 


K; 


-4(n  ;7?Sin  e  +  i  4  (n2)  z2zSin0 


o 


IT 


Thus  our  constraint  becomes 


(Ar ) << 


2  4 (n2-l-Sin2  6 ) 


(31) 


where  we  define  the  "angle"  0  by 


tan  0  = 


K  FFT(d) 

Ko 

FFT  (ij' ) 

(32) 


K  is  the  spacial  Fourier  transform  variable  from  Z  to  K, 
and  the  (n2)zzzz  and  (n2)zzz  terms  have  been  dropped. 

The  depth  mesh  increment  is  constrained  by  requiring 
the  initial  mesh  spacing  to  be  less  than  the  "width"  of 
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the  Gaussian  beam  initial  condition.  As  the  solution  is 


marched  forward,  the  spectrum  is  partitioned  into  4  equal 
bins.  The  energy  Ej  in  each  bin  is  computed  and  the 
following  energy  ratios  are  examined: 


Ei+  E2+  E3+  E4 

e3  +  E4 
E1  ^2 


1.  If  -20dB<Ri<-14dB  a  transform  aliasing  warning  is 
issued.  Five  such  warnings  will  terminate  the  run. 

2.  If  P.2^.-14dB  abiasing  is  considered  severe  and  the 
run  is  immediately  terminated. 

3.  If  R2<-70dB  and  R3<-60dB,  the  field  is  considered 
to  be  oversampled  tnd  the  transform  size  is  reduced. 

Our  final  numerical  algorithm  then  consists  of  the 
following  steps: 

1.  Fourier  sine  transform  ^ 

2.  Compute  the  energy  ratios  and  the  "RMS  angle"  9  . 

3.  Compute  a  ^rnew  at  each  mesh  point  using  the  range 
step  constraint  and  search  for  the  minimum. 

4.  Apply  the  operator 

i(Arold  Arnew)  A  ^ 


16 


5.  Fourier  sine  transform  ij; 

6.  Apply  the  operator 

iArnev  B, 

e  ip 

7.  Advance  the  range  r=rolc2  +  Arnew 
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IV 


USE  AND  LIMITATIONS  OF  THE  MODEL 


Most  of  the  limitations  of  the  parabolic  equation 
underwater  acoustic  propagation  model  are  imposed  by  the 
acoustic  frequency  or  are  associated  with  modeling  the 
ocean  bottom.  The  assumption  that  the  propagation  is  at 
small  angles  to  a  preferred  direction  in  a  medium  with 
slowly  varying  inhomogenities  is  implicit  in  the  para¬ 
bolic  approximation  (equation  7) .  In  addition  other 
limitations  are  imposed  by  the  numerical  algorithm  used 
here.  It  is  evident  from  the  third  order  term  (equation 
30)  that  the  allowable  range  step  is  limited  by  the  sound 
speed  gradients.  Large  gradients  such  as  those  found  in 
bottom  sediment  layers  imply  small  range  steps  and  can 
lead  to  unacceptable  computation  times.  The  computer 
code  will  issue  a  warning  if  the  predicted  range  step  is 
less  than  the  acoustic  wavelength.  Five  such  warnings 
will  cause  termination  of  the  calculation. 

A  small  sound  speed  gradient  approximation  is  also 
made  in  the  solution  of  the  integral  equation  that  leads 
to  the  parabolic  phase  velocity  correction  mapping 
(equation  12) .  Hence  the  mapping  presented  here  is 
applicable  only  to  modes  with  turning  points  in  the  water 


column. 


Since  the  required  depth  mesh  spacing  is  related  to 
the  acoustic  frequency,  a  high  frequency  limit  is  imposed 
by  the  maximum  size  transform  that  can  be  handled  by  the 
computer  code.  The  transform  size  is  in  turn  limited  by 
the  available  storage  and  computation  time  constraints. 
The  time  per  ranqe  step  is  proportional  to  the  transform 
size  and  hence  to  the  frequency.  Higher  frequencies  also 
demand  smaller  range  steps.  Thus,  computation  time  for 
a  given  oceanographic  environment  is  roughly  proportional 
to  the  square  of  the  frequency. 

Fortunately  in  most  ocean  environments,  long  range 
low  frequency  sound  propagation  is  dominated  by  rays  at 
small  grazing  angles.  Energy  at  steeper  angles  is  lost 
due  to  bottom  absorption  and  radiation.  In  these  circum¬ 
stances  the  parabolic  equation  technique  is  suitable  for 
modeling  the  refracted  and  refracted  surface-reflected 
paths  that  account  for  the  propagation  of  energy  to  long 
ranges.  Since  the  parabolic  equation  model  is  a  wave 
model,  diffraction  and  all  other  full  wave  effects  are 
automatically  included.  In  shallow  water,  however,  or 
in  the  deep  ocean  when  bottom  interaction  is  important, 
the  parabolic  equation  model  in  its  present  form  is  not 
appropriate . 
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V  THE  COMPUTER  CODE 

The  PE  computer  program  consists  of  a  main  program 
driver  and  14  computational  subroutines.  In  addition, 
two  environmental  input  subroutines  must  be  supplied  by 
the  user.  Output  is  written  to  an  auxiliary  (tape  or  disk) 
file  as  individual  records  containing  the  range  in  nautical 
miles  and  transmission  loss  (dB  re  1  yard)  at  up  to  2C 
output  depths.  Optional  output  includes  a  printed  trans¬ 
mission  loss  table  and  a  line  printer  field  plot. 

The  user  must  supply  two  subroutines: 

FORTRAN  FUNCTION  ZB(R) 

FORTRAN  SUBROUTINE  SVP (NC , Z , C , RNFXT) 

ZB  returns  the  bottom  depth  in  feet  at  range  R  in  feet. 

SVP  returns  an  NC  point  (£100)  sound  velocity  profile 

table  in  arrays  Z  and  C  and  RNEXT ,  the  range  of  the  next 

profile  on  the  track  in  feet.  Z  is  a  depth  array  in  feet 

or  meters.  C  is  a  sound  speed  array  in  ft/sec  or  m/sec. 

The  core  requirement,  excluding  the  two  user  supplied 
environmental  subroutines,  is  approximately  112,000  octal 
(38,000  decimal)  on  the  CDC  6400/6600/6700. 

The  input  deck  structure  can  be  summarized  as  follows: 

Card  1  FORMAT  (8A10) 

Run  title  (80  columns  text) 

Card  2  FORMAT  (515) 

ND  -  Number  of  output  depths  1<ND<20 
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IFLAT  =  0  Bottom  is  flat 

/  0  Bathymetry  will  he  supplied 
by  user  FUNCTION  ZB(R) 

IPRNT  =0  No  transmission  loss  table 
printed 

/  0  Print  transmission  loss  table 

NPLT  =0  No  line  printer  field  plot 
generated 

_<  120  Generate  line  printer  field 
plot  of  NPLT  depths 

N<6  Program  will  select  transform 

size 

N>6  Program  will  run  at  specified 

transform  size 

The  current  version  of  the  program  is  dimensioned  for  tKll. 
The  maximum  transform  size  is  specified  by  NMAX=2N  in 
SUBROUTINE  PETL.  NMAX  must  correspond  to  the  dimensions 
of  arrays  in  common  blocks  /FIELD/,  /TABLE/,  /WTR/  and 
arrays  B,  JI,  and  ST  in  SUBROUTINE  RST . 

Card  3  FORMAT  (3F10.2) 

ZS  -  Input  depth  in  feet 
F  -  Frequency  in  Hertz 
CQ  -  Reference  sound  speed 

If  C0^o,  the  input  environment  will  be  transformed  to  reduce 
the  parabolic  phase  velocity  error. 

If  C0>o  the  specified  reference  sound  speed  will  be  used  and 

the  parabolic  phase  velocities  will  not  be  corrected. 

Card  4  FORMAT  (7F10.2) 

DMAX  -  The  maximum  water  depth  on  the  track 
in  feet 
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RMAX  -  The  maximum  range  of  the  calculation 
in  nautical  miles 

DR  -  The  range  step  in  nautical  miles 
If  DR<0 ,  the  code  will  select  range  steps  aased  on  the 
split-step  algorithm  truncation  error  estimates.  If 
DR>0 ,  the  code  will  run  at  the  specified  range  step  and 
all  error  checks  are  disabled. 

CDl  -  Minimum  depth  for  the  line  printer 
field  plot  in  feet 

CD2  -  Maximum  depth  for  the  line  printer 
field  plot  in  feet 

CLMIN  -  Minimum  loss  for  the  line  printer 
field  plot  in  dB 

DCL  -  Loss  increment  for  the  line  printer 
field  plot  in  dB 

Card  5  FORMAT  (3F10.2) 

D (I )  -  ND  output  depths  in  feet 
The  formats  for  the  environmental  inputs  (bathymetry  and 
sound  velocity  profiles)  will  be  determined  by  the  user 
supplied  subroutines . 

The  sequence  of  calculations  proceeds  essentially  as 
follows : 

1.  The  input  data  are  read  and  printed  in  the  main 
program  driver. 

2.  SUBROUTINE  PETL  is  called  to  compute  the  transmission 
loss . 

3.  SUBROUTINE  PETL  calls  user  supplied  SUBROUTINE  SVP 
to  obtain  the  initial  sound  velocity  profile.  The  depth 

mesh  is  determined  if  the  transform  size  has  not  been  specified 
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and  relevant  constants  are  computed.  The  transform  size 
is  checked  to  determine  if  array  dimensions  are  exceeded. 
SUBROUTINE  FILTER  is  called  to  evaluate  the  index  of 
refraction  on  the  depth  mesh  by  linear  interpolation  in 
the  sound  speed  profile  table.  If  the  parabolic  phase 
velociiies  are  to  be  corrected  the  transformed  mesh  is 
generated.  A  1-2-1  smoothing  filtei  is  applied  on  the 
index  of  refraction  mesh  to  reduce  transform  aliasing. 
SUBROUTINE  SOURCE  is  called  to  define  the  initial  field 
as  a  Gaussian  beam  with  the  amplitude  and  standard 
deviation  selected  to  match  asymptotically  the  point 
source  solution  of  the  reduced  wave  equation  in  a  homo¬ 
geneous  medium.  The  range  loop  is  then  initiated  and 
SUBROUTINE  STEP  is  called  to  advance  the  solution  one 
range  step. 

4.  If  the  range  step  is  fixed  SUBROUTINE  STEP 
advances  the  solution  the  specified  range  step  using  the 
Tappert-Hardin  split-step  Fourier  algorithm  and  returns 
to  PETL.  If  the  range  step  is  not  specified,  every  fifth 
call  to  STEP  causes  the  sine  transform  to  be  examined  for 
evidence  of  aliasing.  If  the  ratio  of  energy  in  the  last 
1/4  of  the  spectrum  to  the  total  energy  is  greater  than 
-20  dB,  STEP  will  set  the  error  return  flag  and  a  warning 
message  will  be  printed  upon  returning  to  PETL.  If  aliasing 
is  severe  (>p-14dB)  the  run  will  be  terminated.  On  the  other 
hand,  if  the  spectral  energy  distribution  permits,  the 
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transform  size  will  he  reduced.  STEP  then  computes  a 
new  range  step  at  each  mesh  point  using  the  split-step 
algorithm  truncation  error  estimates .  The  mesh  is 
searched  for  the  minimum  step  size.  If  the  minimum  step 
is  less  than  the  acoustic  wavelength,  the  error  flag  is 
set  and  a  warning  message  will  be  printed  upon  returning 
to  PETL.  The  program  will  then  attempt  to  continue  with 
the  current  range  step.  If  the  new  range  step  is  not  small, 
and  the  relative  change  in  step  size  is  greater  than  25%, 

STEP  will  replace  the  current  range  step  with  the  new  step, 
reconstruct  all  stored  tables,  and  advance  the  solution. 

5.  Upon  returning  to  PETL  from  STEP,  the  range  is 
advanced  and  if  required,  a  new  sound  velocity  profile 

is  obtained  and  filtered.  Transmission  loss  at  the  output 
depths  is  interpolated  from  the  field  mesh  and  written  to 
a  tape  or  disk  file.  If  requested  (NPLT>0)  the  line  printer 
field  plot  is  generated.  The  error  return  flag  from 
SUBROUTINE  STEP  is  examined  and  if  set,  appropriate  warning 
messages  are  printed  and  counted.  The  range  loop  will  be 
continued  until  the  maximum  range  is  attained  or  five  warnings 
have  been  issued.  PETL  then  returns  to  the  main  program 
driver. 

6.  If  requested  (IPRNT^O) ,  the  main  program  prints 
the  transmission  loss  table  and  the  run  is  terminated. 

Small  range  step  warning  usually  result  from  dis¬ 
continuities  or  extreme  gradients  in  the  sound  velocity 


profile  or  an  environment  in  which  the  sound  velocity  is 
not  slowly  varying  in  range.  Large  changes  in  the  sound 
velocity  profile  between  range  steps  cannot  be  handled  by 
the  numerical  algorithm.  Transform  aliasing  warninas 
indicate  that  the  field  is  not  adequately  sampled  in  depth. 
The  code  attempts  to  select  the  proper  transform  size; 
hovrever,  since  the  sampling  required  depends  on  the  solution, 
this  is  not  alvav's  possible.  If  the  code  terminates  because 
of  transform  aliasing,  it  is  necessary  to  set  the  next 
larger  transform  size  on  card  1.  The  components  of  the 
computer  code  are  documented  in  Appendix  A. 
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APPENDIX  A 


DESCRIPTION  OE  THE  CODE  COMPONENTS 
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APPENDIX  A 


PROGRAM  PE 

PE  is  the  innut/output  driver  for  the  parabolic 
equation  model.  It  defines  the  FORTRAN  logical  units, 
reads  and  prints  the  input  data,  calls  PETL  to  create  a 
transmission  loss  file  on  FORTRAN  Unit  LT  (dish  or  tape) , 
checks  for  an  error  return  from  PETL,  and  prints  an  output 
transmission  loss  table  if  desired. 

FORTRAN  logical  units  LC,  I.P,  and  I.T  must  be  defined 
by  the  first  throe  executable  statements  of  PF.  I.C  is 
the  input  device  (card  reader) ,  LP  is  the  output  device 
(line  printer) ,  and  LT  is  the  file  output  device. 

CARD  INPUTS  (FORTRAN  Unit  LC) 

Card  1  -  FORMAT  (8A10) 

TITLE  -  Run  title  (text  10  characters/word) 

Card  2  -  FORMAT  (515) 

ND  -  Number  of  output  depths 
IFLAT  -  Flat  bottom  flap 
IPRNT  -  Print  flag 

NPLT  -  Number  of  line  printer  field  plot  depths 
N  -  Transform  size 
Card  3  -  FORMAT  (3F10.2) 

ZS  -  Input  depth  in  feet 
F  -  Frequency  in  Hertz 

CO  -  Reference  sound  speed  in  ft/sec  or  m/sec 
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Card  4  -  FORMAT  (7F10.2) 

DMAX  -  Maximum  depth  in  feet 
RMAX  -  Maximum  range  in  nautical  miles 
DR  -  Range  step  in  nautical  miles 
CD1  -  Minimum  field  plot  depth  in  feet 
CD2  -  Maximum  field  plot  depth  in  feet 
CLMIN  -  Minimum  field  plot  loss  in  dB 
DCL  -  Field  plot  loss  increment  in  dB 
Card  5  -  FORMAT  (8F10.2) 

D  ( 1 )  -  ND  output  depths  in  feet 
COMMON  OUTPUTS 

/UNITS/  LC  -  FORTRAN  input  unit 

LF  -  FORTRAN  output  unit 
LT  -  FORTRAN  file  output  unit 
/HERTZ/  CO  -  Reference  sound  speed 

F  -  Frequency 

/PLT/  TITLE  -  Run  title 

NPI.T  -  Number  of  field  plot  depths 
CLMIN  -  Minimum  loss 
DCL  -  Loss  increment 
CDl  -  Minimum  plot  depth 
CD2  -  Maximum  plot  depth 
/MESH/  DR  -  Range  step 


N  -  Transform  size 


PRINTED  OUTPUT  (FORTRAN  Unit  LP) 

PE  prints  the  input  data,  and  if  IPP.NT/0,  a  range- 
transmission  loss  table  for  the  output  depths. 

SUBROUTINE  PETL 

PETL  sets  the  maximum  transform  size,  establishes  the 
maximum  depth  sample  in  the  transform,  defines  constants, 
sets  up  the  initial  conditions,  controls  the  range  loop, 
and  creates  an  output  transmission  loss  file.  If  the  trans¬ 
form  size  has  not  been  specified  (N<6) ,  PETL  will  select 
a  transform  size  such  that  the  depth  mesh  increment  is 
less  than  the  beam  width  of  the  Gaussian  beam  initial 
condition.  If  the  reference  sound  speed  has  not  been 
specified  (C0£O) ,  PETL  sets  the  phase  velocity  correction 
flag  to  transform  the  environment  to  reduce  the  parabolic 
phase  velocity  error  (MC-1) ,  and  takes  the  minimum  sound 
speed  as  the  reference  sound  speed.  The  dimensions  of 
arrays  in  common  blocks  /FIELD/  and  /TABLE/  must  correspond 
to  the  maximum  transform  size  established  in  PETL. 

CALLING  PROGRAM 
PE 

PARAMETER  INPUTS 

ZS  -  Input  depth  in  feet 
ND  -  Number  of  output  depths 
D  -  Output  depths  in  feet 


DMAX  -  Maximum  depth  in  feet 
RMAX  -  Maximum  range  in  feet 
IFLAT  -  Flat  bottom  flag 


COMMON  INPUTS 
/UNITS/ 

/HERTZ/ 

/PLT/ 

/MESH/ 


LP  -  FORTRAN  output  unit 

LT  -  FORTRAN  file  output  unit 

CO  -  Reference  sound  speed 

H  -  Mesh  increment  in  transform  space 

HK  -  Ratio  of  mesh  increment  in  transform 
space  to  reference  wave  number 

FK  -  Reference  wave  number 

WL  -  Acoustic  wavelength 

DCD  -  Field  plot  depth  increment 

CD  -  Field  plot  depth  array 

R  -  Current  range  in  feet 

NR  -  Range  step  count 

KR  -  Flagged  step  count 

DZ  -  Depth  mesh  increment  in  feet 

ZMAX  -  Maximum  depth  sample  in  the 
transform  in  feet 

IB  -  Depth  index  of  bottom  interface 
N  -  Transform  size 

NPTS  -  Number  of  field  mesh  points  (2N-1) 
N2  -  1/2  the  number  of  field  mesh  points 
N 4  -  1/4  the  number  of  field  mesh  points 
NL4  -  3/4  the  number  of  field  mesh  points 
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NA  -  Number  of  points  in  the  artificial 
attenuation  table 

\'W  -  Number  of  points  in  the  water  column 

’W  -  Current  bottor  depth  in  feet 

HALF  -  Half  the  number  of  field  mesh  points 

/PHASE/  NC  -  Number  of  points  in  current  sound 

speed  profile 

1  -  Sound  speed  depth  array 

C  -  Sound  speed  array 

MC  -  Phase  velocity  correction  flag 

DM  -  Transformed  output  depths 

/OUTBUF/  RNM  -  Current  range  in  nautical  miles 

TL  -  Transmission  loss  at  the  output 
depths  in  dB 

FILE  OUTPUT 

Unit  LT  -  One  record  for  each  range  step  including  the 
range  in  nautical  miles  and  the  transmission 
loss  at  up  to  20  output  depths 

PRINTED  OUTPUT 

PETL  prints  the  transform  size  and  appropriate  warnings 
if  the  error  return  flag  from  STEP  is  set. 

FUNCTION  TLOSS 

TLOSS  interpolates  the  absolute  value  of  the  field  at 
depth  Z  from  the  field  mesh  and  returns  the  transmission 
loss.  The  interpolation  is  linear  in  the  absolute  value 


of  the  field.  The  maximum  loss  permitted  is  180  dB. 

CALLING  PROGRAMS 
PETL 
FLD 

PARAMETER  INPUTS 

RR  -  Reciprocal  range  in  inverse  feet 
Z  -  Depth  in  feet 
COMMON  INPUTS 

/FIELD/  PR  -  Real  part  of  the  field  mesh 

PI  -  Imaginary  part  of  the  field  mesh 
/MESH/  DZ  -  Depth  mesh  increment  in  feet 

FUNCTION  OUTPUT 

TLOSS  -  Transmission  loss  at  depth  Z 

SUBROUTINE  SOURCE 

SOURCE  constructs  the  field  to  be  used  as  an  initial 
condition  for  the  numerical  integration.  The  field  is 
defined  as  a  Gaussian  beam  at  zero  range.  The  parameters 
of  the  beam  are  selected  to  match  asympotically  the  point 
source  solution  of  the  elliptic  wave  equation.  The  width, 
GW,  is  given  by 

2 

GW  =  — 

FK 

where  FK  is  the  average  wave  number.  The  amplitude,  GA,  is 
given  by 


The  real  part  of  the  field,  PR,  is  set  equal  to  the  beam 
minus  its  image. 


where  ZM  is  the  depth  mesh  point  and  the  exponent  is  con¬ 
strained  to  be  greater  than  -42.  The  imaginary  part  of  the 
field  is  set  to  loro. 

PI  =  0 

CALLING  PROGRAMS 
PETL 

PARAMETER  INPUTS 

ZS  -  Input  depth  in  feet 
COMMON  INPUTS 

/HERTZ/  EK  -  Reference  wave  number 

/MESH/  DZ  -  Depth  mesh  increment  in  feet 

COMMON  OUTPUTS 

/FIELD/  PR  -  Real  part  of  the  initial  field 

PI  -  Imaginary  part  of  the  initial  field 
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FUNCTION  SPEED 


SPEED  linearly  interpolates  the  sound  speed  at  depth 
D 

CALLING  PROGRAMS 
PETL 
FILTER 

PARAMETER  INPUT 

D  -  Depth  in  feet 
COMMON  INPUTS 

/PHASE/  NC  -  Number  of  sound  velocity  profile 

points 

7,  -  Depth  array 
C  -  Sound  speed  array 

FUNCTION  OUTPUT 

SPEED  -  Sound  speed  at  depth  D 

SUBROUTINE  FILTER 

FILTER  evaluates  the  index  of  refraction  on  the  field 
mesh  and  transforms  the  environment  to  reduce  the  phase 
velocity  error  inherent  in  the  parabolic  approximation. 
The  input  sound  speed  profile  is  checked  and  if  required 
the  units  are  converted  to  feet  and  ft/sec.  The  sound 
speed  is  assumed  to  be  piecewise  linear  in  depth.  The 
output  array  FN  contains  n^-l  on  a  uniform  mesh  where  n 
is  the  index  of  refraction  or  the  transformed  index  of 
refraction 
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1/2 

n=  (2n-l) 

if  the  parabolic  phase  velocity  correction  is  applied. 
The  output  array  is  filtered  to  smooth  the  gradient  dis¬ 
continuities  introduced  by  the  boundaries  in  the  input 
profile.  The  transformed  depths  are  giver,  by 

’  co  ]  1/2 

2  =  2  - 

C  (  Z  ) 

The  inverse  mapping 


2C0 


where  g  is  the  sound  speed  gradient 
Zk  is  the  boundary  depth 

C k  is  the  sound  speed  at  the  boundary  depth 
is  used  to  evaluate  the  sound  speed  on  the  transformed  depth 
mesh.  The  flag  MC  is  used  to  signal  the  pahse  velocity 
correction  mapping. 

MC  =  1  Phase  velocity  correction 
MC  =  2  No  phase  velocity  correction 
CALLING  PROGRAM? 

PETL 

PARAMETER  INPUTS 

ND  -  Number  of  output  depths 
D  -  Output  depth  array 
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COMMON  INPUTS 


/UNITS/ 

/HERTZ/ 

/PHASE/ 


/PLT 


/MESH/ 


LP  -  FORTRAN  output  unit 

CO  -  Reference  sound  speed  in  ft/sec 

NC  -  Number  of  sound  speed  profile  points 

Z  -  Depth  array 

C  -  Sound  speed  array 

MC  -  Phase  velocity  correction  flag 

NPLT  -  Number  of  field  plot  depths 

CDl  -  Minimum  field  plot  depth 

DCD  -  Field  plot  depth  increment 

DZ  -  Depth  mesh  increment  in  feet 

MW  -  Maximum  number  of  points  in  the  water 
column 

NPTS  -  Number  of  field  mesh  points 


COMMON  OUTPUTS 

/PHASE/  DM  -  Transformed  output  depths 

/PLT/  CD  -  Transformed  field  plot  depths 

/TABLE/  FN  -  Filtered  index  of  refraction  array 

PRINTED  OUTPUT 

Warning  message  if  the  sound  speed  profile  is  not  defined 
to  the  surface.  In  this  case  the  run  is  aborted. 


SUBROUTINE  SET 

SET  constructs  all  tables  that  are  a  function  of  the 
range  step.  The  artificial  bottom  attenuation  table  is  of 
the  form 


36 


-0 . 01*DR*e 


I 

I 


DZ*I-DZ*NA\ 

173  f Z*NA / 


where  DZ  is  the  depth  increment  in  feet 
DR  is  the  range  step  in  feet 
NA  is  the  number  of  points  in  the  table 
The  second  derivative  transform  table  is  of  the  form 


S 


DR 


2Fk 


K' 


e 

2N-1 


where  DR  is  the  range  step  in  feet 
FK  is  the  average  wave  number 

K  is  the  spacial  Fourier  transform  variable  from 
Z  to  K 

N  is  the  transform  size 

The  factor  2N_-1-  provides  the  transform  normalization.  Th< 
real  part  is  returned  in  array  SR.  The  imaginary  part  is 
returned  in  array  SI. 

SUBROUTINE  INDEX  is  called  to  generate  an  index  of 
refraction  table  of  the  form 
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where  DR  is  the  range  step  in  feet 
FK  is  the  average  wave  number 
n  is  the  index  of  refraction 

The  real  part  is  returned  in  array  UR.  The  imaginary  part 
is  returned  in  array  UI . 

CALLING  PROGRAMS 
PETL 
STEP 

COMMON  INPUTS 

/HERTZ/  H  -  Mesh  increment  in  transform  space 

FK  -  Reference  wave  number 

/MESH/  DR  -  Range  step  in  feet 

DZ  -  Depth  increment  in  feet 

NPTS  -  Number  of  depth  mesh  points 

NA  -  Number  of  points  in  the  attenuation 
table 

HALF  -  2n_1  where  N  is  the  transform  size 

COMMON  OUTPUTS 

/TABLE/  A  -  Attenuation  table 

SR  -  Real  part  of  the  transform  table 

SI  -  Imaginary  part  of  the  transform  table 

UR  -  Real  part  of  the  index  of  refraction 
table 
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/HERTZ/ 


UI  -  Imaginary  part  of  the  index  of 
refraction  table 

FACTOR  -  Constant  for  generating  index 
of  refraction  table  (DR*FK/2) 


SUBROUTINE  INDEX 

INDEX  constructs  an  index  of  refraction  trable  of  the 


form 


DR*FK  ,  2 
2  (n  1} 


U  =  e 


where  DR  is  the  range  step 

FK  is  the  reference  wave  number 
n  is  the  index  of  refraction 

The  real  part  of  the  table  is  returned  in  array  UR.  The 
imaginary  part  is  returned  in  array  UI.  The  artificial 
attenuation  is  incorporated  in  the  output  table 
CALLING  PROGRAMS 
PETL 


COMMON  INPUTS 


/HERTZ/ 


/MESH/ 


FACTOR  -  DR*FK/2  where  DR  is  the  range  step 
and  FK  is  the  reference  wave  number 

NPTS  -  Number  of  depth  mesh  points 

NW  -  Maximum  number  of  points  in  the  water 
column 

NA  -  Number  of  points  in  the  attenuation 
table 

IB  -  Mesh  index  for  start  of  bottom 


39 


/TABLE/  A  -  Attenuation  table 

COMMON  OUTPUTS 

/TABLE/  UR  -  Real  part  of  index  of  refraction 

table 

UI  -  Imaginary  part  of  index  of  refraction 
table 

SUBROUTINE  STEP 

STEP  advances  the  field  one  range  step  using  the 
Tappert-Hardin  split-step  Fourier  algorithm.  If  the  range 
step  has  not  been  fixed,  STEP  will  perform  error  checks 
and  use  the  numerical  algorithm  trunfation  error  estimates 
to  compute  a  new  range  step  every  fifth  call  (i.e.  when 
NR=KR) .  If  the  relative  change  in  step  size  is  greater 
than  25%,  the  current  range  step  is  updated  and  new  tables 
are  constructed.  If  the  predicted  range  step  is  less  than 
the  acoustic  wavelength  or  the  ratio  of  the  energy  in  the 
last  1/4  of  the  spectrum  to  the  total  energy  is  greater 
than  -20  dB  (evidence  of  transform  aliasing) ,  the  error 
return  flag  is  set  (FLAG  =  0) .  In  addition  two  other 
spectral  energy  ratios  are  examined  and  if  all  criteria  are 
met,  the  field  is  considered  to  be  oversampled  in  depth  and 
the  transform  size  is  reduced. 

CALLING  PROGRAMS 
PETL 

COMMON  INPUTS 

/FIELD/  PR  -  Real  part  of  field  at  range  R 

PI  -  Imaginary  part  of  field  at  range  R 


/HERTZ/ 


/MESH/ 


/TABLE/ 


II  -  Mesh  increment  in  transform  space 

HK  -  Ratio  of  mesh  increment  in  transform 
space  to  the  reference  wave  number 

FK  -  Reference  wave  number 

WL  -  Acoustic  wavelength 

DR  -  Current  range  step  in  feet 

UR  -  Current  range  step  count 

KR  -  Current  flagged  step  count 

DZ  -  Depth  mesh  increment  in  feet 

IB  -  Mesh  index  for  start  of  Dottom 

N  -  Transform  size 

NPTS  -  Number  of  field  mesh  points 

H2  -  1/2  the  number  of  field  mesh  points 

U4  -  1/4  the  number  of  field  mesh  points 

NL4  -  3/4  the  number  of  field  mesh  points 

NA  -  Number  of  attenuation  table  points 

NW  -  Number  of  mesh  points  in  the  water 
column 

HALF  -  1/2  the  number  of  field  mesh  points 

A  -  Artificial  attenuation  table 

FN  -  n^-1  where  n  is  the  index  of  refraction 

SR  -  Real  part  of  transform  table 

SI  -  Imaginary  part  of  transform  table 

UR  -  Real  part  of  index  of  refraction  table 

UI  -  Imaginary  part  of  index  of  refraction 
table 
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PARAMETER  OUTPUTS 


FLAG  -  Error  return  flag 

COMMON  OUTPUTS 

/FIELD/  PR  -  Real  part  of  field  at  range  R+DR 

PI  -  Imaginary  part  of  field  at  range 
R+DR 

/MESH/  DR  -  Updated  range  step 

KR  -  Updated  flagged  step  count 
DZ  -  Depth  mesh  increment 
N  -  Transform  size 

NPTS  -  Number  of  field  mesh  points 

N2  -  1/2  the  number  of  field  mesh  points 

N4  -  1/4  the  number  of  field  mesh  points 

NL4  -  3/4  the  number  of  field  mesh  points 

NA  -  Number  of  points  in  the  attenuation 
table 

NW  -  Maximum  number  of  mesh  points  in  the 
water  column 

HALF  -  1/2  the  number  of  field  mesh  points 

/TABLE/  A  -  Artificial  attenuation  table 

FN  -  n^-i  where  n  is  the  index  of  refraction 

SR  -  Real  part  of  the  transform  table 

SI  -  Imaginary  part  of  the  transform  table 

UR  -  Real  part  of  the  index  of  refraction 
table 

UI  -  Imaginary  part  of  the  index  of  refraction 
table 
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FUNCTION  GAUSS 


GAUSS  evaluates  the  Gaussian  function 


GA  *  e 


The  exponent  is  not  allowed  to  be  less  than  -42. 
CALLING  PROGRAMS 
SET 
SOURCE 

PARAMETER  INPUTS 

GA  -  Gaussian  amplitude 
GD  -  Gaussian  depth  in  feet 
GW  -  Gaussian  width  in  feet 
Z  -  Depth  in  feet 

FUNCTION  OUTPUT 


GAUSS 


SUBROUTINE  FLD 

FLD  is  a  routine  for  plotting  the  transmission  loss 
field  generated  by  the  parabolic  equation  model  on  the 
line  printer.  Symbols  corresponding  to  five  transmission 
loss  bins  at  up  to  120  depths  are  plotted  at  approximately 
one  nautical  mile  range  increments.  If  the  plot  depth  is 
greater  than  the  current  water  depth,  B  is  printed.  A 
transmission  loss  scale  and  a  depth  scale  are  printed 
the  first  time  FLD  is  called. 


CALLING  PROGRAMS 


PETL 

PARAMETER  INPUTS 


COMMON  INPUTS 
/UNITS/ 
/OUTBUF/ 
/PLT/ 


/MESH/ 


RR  -  Reciprocal  range 

LP  -  FORTRAN  output  unit 

RNM  -  Current  range  in  nautical  miles 

TITLE  -  Run  title  (text) 

NPLT  -  Number  of  plot  depths 

LCR  -  Last  range  plotted 

CLMIN  -  Minimum  plot  loss  in  dB 

DCL  -  Loss  increment  in  dB 

CD1  -  Minimum  plot  depth  in  feet 

DCD  -  Plot  depth  increment  in  feet 

CD  -  Plot  depth  array  in  feet 

R  -  Current  range  in  feet 

DR  -  Current  range  step  in  feet 

ZW  -  Current  water  depth  in  feet 


COMMON  OUTPUTS 

/PLT/  LCR  -  Current  plot  range 

PRINTED  OUTPUT 

Line  printer  plot  of  current  range  and  transmission 


loss  bins  at  up  to  120  depths. 


SUBROUTINE  RST 


RST  is  a  sine  transform  driver  for  a  real  vector 
mixed  radix  (8-4-2)  fast  Fourier  synthesis  algorithm 
developed  by  Bergland.^  The  sine  transform  algorithm  is 
from  Cooley,  Lewis  and  Welsh. ®  RST  defines  constrants , 
sets  up  the  necessary  tables,  forms  the  input  Fourier 
coefficients  for  the  fast  Fourier  synthesis  subroutines, 
and  forms  the  sine  transform  from  the  results.  Upon 
returning  from  RST  the  real  input  vector  has  been  replaced 
by  its  finite  discrete  sine  transform.  The  dimensions  of 
arrays  in  RST  must  correspond  to  the  largest  vector  to  be 
transformed . 

ARRAY  REQUIRED  DIMENSIONS 


B 

ST 

JI 

CS 

SS 


2n 

2n-l 

2n~l  -1 

2n-4 
2n~ 4  _i 


where  the  real  input  vector  is  of  dimension  2n-l 


CALLING  PROGRAMS 
STEP 

PARAMETER  INPUTS 

X  -  Real  vector 

N  -  Transform  size  (input  vector  of  size 
2N-1) 


PARAMETER  OUTPUTS 


X  -  Sine  transform  of  the  input  vector 


COMMON  OUTPUTS 


Size  of  trigonometric  table  for 
a  given  radix  8  iteration 

Cosine  table  stored  in  bit-reversed 
order 

SS  -  Sine  table  stored  in  bit-reversed 
order 

SUBROUTINE  R8SYN 

R8SYN  performs  one  pass  of  a  radix  8  real  vector  fast 
Fourier  synthesis. 

CALLING  PROGRAM 
RST 

PARAMETER  INPUTS 

INT  -  Length  of  current  pass 

BO 

Bl 

B2 

B3  -  Input  vectors  for  current  radix 
8  pass 
B4 

B5 

B6 

B7 

COMMON  INPUTS 

/WTS/  NT  -  Length  of  trigonometric  table 

CS  -  Cosine  table  stored  in  bit-reversed 
order 

SS  -  Sine  table  stored  in  bit-reversed  order 


/WTS/  NT  - 

CS  - 


PARAMETER  OUTPUTS 


BO 

Bl 

B2 

B3 

B4  -  Output  vectors  of  current  radix 
8  pass 
B5 

B6 

B7 


SUBROUTINE  R4SYN 

R4SYN  performs  a  single  pass  radix  4  real  vector  fast 
Fourier  synthesis. 

CALLING  PROGRAMS 
RST 

PARAMETER  INPUTS 

INT  -  Length  of  radix  4  pass 
BO, 

B! 

B2 l-  Input  vectors  for  radix  4  pass 
B3  ' 

PARAMETER  OUTPUTS 

BO 

B1 

B2\-  Output  vectors  from  radix  4  pass 
B3  ' 
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SUBROUTINE  R2TR 


R2TR  performs  a  single  pass  radix  2  fast  Fourier 
transform. 

CALLING  PROGRAM 
RST 

PARAMETER  INPUTS 

INT  -  Length  of  radix  2  pass 
BO 

-  Input  vectors  for  radix  2  pass 

'  B1 

PARAMETER  OUTPUTS 

BO 

-  Output  vectors  from  radix  2  pass 
Bl 
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